Administrative boundary regionizer
import geopandas as gpd
from osmnx import geocode_to_gdf
import plotly.express as px
from shapely.geometry import Point, box
from srai.regionizers import AdministrativeBoundaryRegionizer
/opt/hostedtoolcache/Python/3.8.14/x64/lib/python3.8/site-packages/spherical_geometry/great_circle_arc.py:365: RuntimeWarning: invalid value encountered in divide return P / l
Regionize city¶
Basic usage of the AdministrativeBoundaryRegionizer with a city boundary.
Here admin_level equal to 9 defines city districts in Poland.
wroclaw_gdf = geocode_to_gdf(query=["R451516"], by_osmid=True)
wroclaw_gdf.plot()
<AxesSubplot: >
abr = AdministrativeBoundaryRegionizer(admin_level=9)
wro_result_gdf = abr.transform(gdf=wroclaw_gdf)
wro_result_gdf.head()
Loading boundaries: 0it [00:00, ?it/s][overpass] downloading data: [timeout:60][out:json];(relation["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="9"](51.0426686,16.8073393,51.2100604,17.1762192);); out ids tags; Loading boundaries: 48it [01:04, 1.33s/it]
| geometry | |
|---|---|
| region_id | |
| Osiedle Leśnica | POLYGON ((16.90228 51.18284, 16.90247 51.18202... |
| Osiedle Jerzmanowo-Jarnołtów-Strachowice-Osiniec | POLYGON ((16.83783 51.12350, 16.83855 51.12365... |
| Osiedle Pracze Odrzańskie | POLYGON ((16.91883 51.19454, 16.92015 51.19200... |
| Osiedle Żerniki | POLYGON ((16.93916 51.12235, 16.93851 51.12135... |
| Osiedle Maślice | POLYGON ((16.91968 51.15787, 16.92075 51.15993... |
fig = px.choropleth_mapbox(
wro_result_gdf,
geojson=wro_result_gdf,
color=wro_result_gdf.index,
locations=wro_result_gdf.index,
center={"lat": 51.125, "lon": 16.99},
mapbox_style="carto-positron",
zoom=10.5,
)
fig.update_layout(margin={"r": 0, "t": 0, "l": 0, "b": 0})
fig.update_traces(marker={"opacity": 0.6}, selector=dict(type="choroplethmapbox"))
fig.update_traces(showlegend=False)
minx, miny, maxx, maxy = wroclaw_gdf.geometry[0].bounds
fig.update_geos(
projection_type="equirectangular",
lataxis_range=[miny - 0.1, maxy + 0.1],
lonaxis_range=[minx - 0.1, maxx + 0.1],
showlakes=False,
showcountries=False,
showframe=False,
resolution=50,
)
fig.update_layout(height=600, width=800, margin={"r": 0, "t": 0, "l": 0, "b": 0})
fig.show(renderer="png") # replace with fig.show() to allow interactivity
Regionize country¶
How to return an empty region covering water bodies outside of the land.
Here admin_level equal to 4 defines country regions in Madagascar.
madagascar_bbox = box(
minx=43.2541870461, miny=-25.6014344215, maxx=50.4765368996, maxy=-12.0405567359
)
madagascar_bbox_gdf = gpd.GeoDataFrame({"geometry": [madagascar_bbox]}, crs="EPSG:4326")
abr = AdministrativeBoundaryRegionizer(admin_level=4, return_empty_region=True)
madagascar_result_gdf = abr.transform(gdf=madagascar_bbox_gdf)
madagascar_result_gdf.tail()
Loading boundaries: 0it [00:00, ?it/s][overpass] downloading data: [timeout:60][out:json];(relation["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="4"](-25.6014344215,43.2541870461,-12.0405567359,50.4765368996);); out ids tags; Loading boundaries: 25it [00:45, 1.83s/it]
| geometry | |
|---|---|
| region_id | |
| Atsinanana | POLYGON ((47.67041 -20.06128, 47.67343 -20.044... |
| Diana Region | MULTIPOLYGON (((48.35982 -13.60568, 48.35872 -... |
| Analanjirofo | MULTIPOLYGON (((48.90954 -16.67829, 48.92118 -... |
| Sava | POLYGON ((49.06377 -14.19881, 49.06448 -14.198... |
| EMPTY | MULTIPOLYGON (((50.47654 -25.60143, 45.21301 -... |
fig = px.choropleth(
madagascar_result_gdf,
geojson=madagascar_result_gdf.geometry,
locations=madagascar_result_gdf.index,
color=madagascar_result_gdf.index,
color_continuous_scale=px.colors.sequential.Viridis,
)
fig.update_traces(marker={"opacity": 0.6}, selector=dict(type="choropleth"))
fig.update_layout(coloraxis_showscale=False)
fig.update_traces(showlegend=False)
minx, miny, maxx, maxy = madagascar_bbox.bounds
fig.update_geos(
projection_type="equirectangular",
lataxis_range=[miny - 0.1, maxy + 0.1],
lonaxis_range=[minx - 0.1, maxx + 0.1],
visible=False,
showlakes=False,
showcountries=False,
showframe=False,
resolution=50,
)
fig.update_layout(height=800, width=420, margin={"r": 0, "t": 0, "l": 0, "b": 0})
fig.show(renderer="png") # replace with fig.show() to allow interactivity
Regionize Europe¶
Option to slightly increase the value of toposiplify to simplify geometries even more.
Here admin_level equal to 2 defines countries.
eu_bbox = box(minx=-10.478556, miny=34.633284672291, maxx=34.597916, maxy=70.096054)
eu_bbox_gdf = gpd.GeoDataFrame({"geometry": [eu_bbox]}, crs="EPSG:4326")
abr = AdministrativeBoundaryRegionizer(admin_level=2, toposimplify=0.0005)
eu_result_gdf = abr.transform(gdf=eu_bbox_gdf)
Loading boundaries: 0it [00:00, ?it/s][overpass] downloading data: [timeout:60][out:json];(relation["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="2"](34.633284672291,-10.478556,70.096054,34.597916);); out ids tags; Loading boundaries: 55it [01:52, 2.05s/it]
eu_result_gdf.head()
| geometry | |
|---|---|
| region_id | |
| Netherlands | MULTIPOLYGON (((3.08610 51.57045, 3.09529 51.5... |
| Portugal | POLYGON ((-9.05070 41.86519, -8.87097 41.86393... |
| Morocco | POLYGON ((-5.88545 35.96488, -5.66652 35.93567... |
| Ireland | POLYGON ((-6.71943 55.28017, -6.82172 55.21647... |
| Faroe Islands | POLYGON ((-8.11704 62.11484, -8.11273 62.13026... |
fig = px.choropleth(
eu_result_gdf,
geojson=eu_result_gdf.geometry,
locations=eu_result_gdf.index,
color=eu_result_gdf.index,
color_continuous_scale=px.colors.sequential.Viridis,
)
fig.update_traces(marker={"opacity": 0.6}, selector=dict(type="choropleth"))
fig.update_layout(coloraxis_showscale=False)
fig.update_traces(showlegend=False)
minx, miny, maxx, maxy = eu_bbox.bounds
fig.update_geos(
projection_type="equirectangular",
lataxis_range=[miny - 1, maxy + 1],
lonaxis_range=[minx - 1, maxx + 1],
showlakes=False,
showcountries=False,
showframe=False,
resolution=50,
)
fig.update_layout(height=800, width=1000, margin={"r": 0, "t": 0, "l": 0, "b": 0})
fig.show(renderer="png") # replace with fig.show() to allow interactivity
Toposimplify differences¶
Shows differences in simplification of small regions using four values: 1e-4, 1e-3, 1e-2 and 0.1. Those values are in degress, since it uses Douglas-Peucker simplification algorithm.
1e-4 is the default value and is equal to about 11.1m accuracy.
More info: https://github.com/mattijn/topojson
Here admin_level equal to 6 defines city districts in Singapore.
singapore_bbox = box(minx=103.5111238, miny=1.1263707, maxx=104.1313374, maxy=1.4787511)
singapore_bbox_gdf = gpd.GeoDataFrame({"geometry": [singapore_bbox]}, crs="EPSG:4326")
results = {}
for value in [0.0001, 0.001, 0.01, 0.1]:
abr = AdministrativeBoundaryRegionizer(admin_level=6, toposimplify=value)
results[value] = abr.transform(gdf=singapore_bbox_gdf)
Loading boundaries: 0it [00:00, ?it/s][overpass] downloading data: [timeout:60][out:json];(relation["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="6"](1.1263707,103.5111238,1.4787511,104.1313374);); out ids tags; Loading boundaries: 8it [00:12, 1.51s/it] Loading boundaries: 8it [00:00, 101.83it/s] Loading boundaries: 8it [00:00, 102.11it/s] Loading boundaries: 8it [00:00, 104.03it/s]
minx, miny, maxx, maxy = singapore_bbox.bounds
for epsilon, result in results.items():
fig = px.choropleth_mapbox(
result,
geojson=result,
color=result.index,
locations=result.index,
center={"lat": 1.3119350704252704, "lon": 103.82412242562575},
mapbox_style="carto-positron",
zoom=9.5,
)
fig.update_layout(margin={"r": 0, "t": 0, "l": 0, "b": 0})
fig.update_traces(marker={"opacity": 0.6}, selector=dict(type="choroplethmapbox"))
fig.update_traces(showlegend=False)
fig.update_geos(
projection_type="equirectangular",
lataxis_range=[miny - 0.1, maxy + 0.1],
lonaxis_range=[minx - 0.1, maxx + 0.1],
showlakes=False,
showcountries=False,
showframe=False,
resolution=50,
)
size = len(result.to_json().encode("utf-8"))
fig.update_layout(
height=450,
width=700,
margin={"r": 0, "t": 50, "l": 0, "b": 0},
title_text=f"Toposimplify value: {epsilon} ({size/1000} KB)",
)
fig.show(renderer="png") # replace with fig.show() to allow interactivity
Regionize points¶
How to return original regions without clipping and select them using list of points. Showed using list of metro stations in Paris.
Here admin_level equal to 8 defines communes in France.
import requests
r = requests.get("https://raw.githubusercontent.com/w8r/paris-metro-graph/master/metro.json").json()
stations_gdf = gpd.GeoDataFrame(
{"geometry": [Point(s["longitude"], s["latitude"]) for s in r["nodes"]]}, crs="EPSG:4326"
)
stations_gdf
| geometry | |
|---|---|
| 0 | POINT (2.30159 48.81081) |
| 1 | POINT (2.29703 48.81534) |
| 2 | POINT (2.24947 48.88824) |
| 3 | POINT (2.27210 48.88106) |
| 4 | POINT (2.28937 48.87558) |
| ... | ... |
| 309 | POINT (2.33812 48.88228) |
| 310 | POINT (2.34955 48.88729) |
| 311 | POINT (2.35482 48.84599) |
| 312 | POINT (2.29583 48.87493) |
| 313 | POINT (2.40655 48.87706) |
314 rows × 1 columns
abr = AdministrativeBoundaryRegionizer(admin_level=8, return_empty_region=False, clip_regions=False)
paris_districts_result = abr.transform(gdf=stations_gdf)
Loading boundaries: 0it [00:00, ?it/s][overpass] downloading data: [timeout:60][out:json];is_in(48.8108097,2.3015874);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 1it [00:02, 2.58s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8153367,2.2970255);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 2it [00:04, 2.41s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8882396,2.2494747);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 3it [00:07, 2.51s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8810615,2.2720961);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 4it [00:11, 3.15s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8755816,2.2893655);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 5it [00:14, 3.21s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8336468,2.2435931);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 6it [00:18, 3.48s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8241474,2.2731713);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 7it [00:21, 3.26s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8445799,2.438912);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 8it [00:24, 3.02s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.846435,2.4184422);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 9it [00:26, 2.82s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8653357,2.416682);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 10it [00:29, 2.87s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.9163969,2.2946541);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 11it [00:31, 2.64s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8920499,2.2850318);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 12it [00:33, 2.52s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.857946,2.4357587);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 13it [00:36, 2.45s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.9032278,2.3059879);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 14it [00:38, 2.51s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8264064,2.406201);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 15it [00:41, 2.51s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.89128,2.4025383);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 16it [00:44, 2.80s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8797602,2.4164717);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 17it [00:47, 2.72s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.906204,2.331871);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 18it [00:50, 2.73s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.9231172,2.2862019);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 19it [00:52, 2.74s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.9142573,2.4038072);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 20it [00:56, 2.91s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.9207587,2.4107677);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 21it [00:59, 2.99s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.7963514,2.449235);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 22it [01:01, 2.79s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8033311,2.4456085);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 23it [01:07, 3.68s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.9195753,2.3433012);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 24it [01:10, 3.61s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8051317,2.3637903);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 25it [01:14, 3.44s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.810328,2.3622359);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 26it [01:16, 3.15s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.8159011,2.3773455);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 27it [01:18, 2.94s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.9064052,2.4491919);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 28it [01:22, 3.06s/it][overpass] downloading data: [timeout:60][out:json];is_in(48.818234,2.3195203);area._["boundary"="administrative"]["type"~"boundary|multipolygon"]["admin_level"="8"];relation(pivot); out ids tags; Loading boundaries: 29it [01:25, 2.96s/it]
paris_districts_result.head()
| geometry | |
|---|---|
| region_id | |
| Châtillon | POLYGON ((2.27988 48.81091, 2.28325 48.81056, ... |
| Malakoff | POLYGON ((2.30132 48.82513, 2.31413 48.82226, ... |
| Puteaux | POLYGON ((2.25372 48.88691, 2.25244 48.88551, ... |
| Neuilly-sur-Seine | POLYGON ((2.24562 48.87636, 2.24777 48.87849, ... |
| Paris | POLYGON ((2.23208 48.86951, 2.24046 48.87189, ... |
# Paris
fig = px.choropleth_mapbox(
paris_districts_result,
geojson=paris_districts_result,
color=paris_districts_result.index,
locations=paris_districts_result.index,
center={"lat": 48.858, "lon": 2.353},
mapbox_style="carto-positron",
zoom=10.8,
)
fig2 = px.scatter_mapbox(stations_gdf, lat=stations_gdf.geometry.y, lon=stations_gdf.geometry.x)
fig.add_trace(fig2.data[0])
fig.update_layout(margin={"r": 0, "t": 0, "l": 0, "b": 0})
fig.update_traces(showlegend=False)
fig.update_traces(marker={"opacity": 0.6}, selector=dict(type="choroplethmapbox"))
fig.update_traces(marker_color="black", marker_size=5, selector=dict(type="scattermapbox"))
fig.update_layout(height=800, width=800, margin={"r": 0, "t": 0, "l": 0, "b": 0})
fig.show(renderer="png") # replace with fig.show() to allow interactivity